clc;
clear;
%% 场景建立
startTime = datetime(2023,11,7,8,30,0);
%startTime = datetime(2023,10,29,0,5,34);
stopTime = startTime + seconds(3600*1);
sampleTime = 60;  %60s
timeStep = seconds(sampleTime);                    % 时间步长
timeRange = startTime : timeStep : stopTime;    % 时间轴

sc = satelliteScenario(startTime,stopTime,sampleTime);
%% 地面测站
groundLat = [39.908,46.800,39.505,25.027,18.313,24.480,34.445];  
groundLon = [116.420,130.320,75.929,102.796,109.311,118.09,109.494];
groundAlt = 1000.*[0.049,0.101,1.255,1.955,0.018,0.015,0.557];
numStation = length(groundLon);
gSation = groundStation(sc,'Latitude',groundLat,'Longitude',...
    groundLon,'Altitude',groundAlt);
%% 北斗卫星
satfile = 'starlink_11_7.tle';
sat1 = satellite(sc,satfile);
time_size = length(sat.timerange);
my_sat = sat.sat(1);
%% 星下点轨迹
satpos = zeros(time_size,3);
for i = 1:time_size
    [pos,~] = states(my_sat,sat.timerange(i),"CoordinateFrame","geographic");
    satpos(i,:) = reshape(pos,3,1)';
end
lat = satpos(:,1);
lon = satpos(:,2);
numLon = length(lon);
for i = 1:numLon-1
    if abs(lon(i)-lon(i+1)) >100   %检测经度跳变
        lat = [lat(1:i);NaN;lat(i+1:end)];
        lon = [lon(1:i);NaN;lon(i+1:end)];
    end
end
geoplot(lat, lon,'LineWidth',2,'Color','yellow');
hold on
geoscatter(groundLat,groundLon,"m","d")
geolimits([-90 90],[-180 180]);
geobasemap topographic
title('星下点轨迹');
